R Notebook

if(!require("plotly")) {install.packages("plotly")}

# install.packages("latex2exp")
# install.packages("BiocManager") 
# install.packages("corrplot")
# BiocManager::install("EBImage")

if(!require("lme4")){install.packages("lme4")}
if(!require("lmerTest")){install.packages("lmerTest")}
if(!require("nlme")){install.packages("nlme")}
if(!require("formattable")){install.packages("formattable")}
if(!require("xgboost")){install.packages("xgboost")}
if(!require("processx")) {install.packages("processx")}

library(plotly)
library(lme4)
library(lmerTest)
library(nlme)
library(formattable)
library(xgboost)

### Load libraries
library(EBImage)
library(ggplot2)
library(stringr)
library(gridExtra)
library(latex2exp)
packageVersion('plotly')

[1] ‘4.9.1’

Sys.setenv("plotly_username"="thuynh32")
Sys.setenv("plotly_api_key"="xcSv1yzujDc1IGEwQlr2")

colorBlue = "#007fff"
colorRed = "#ff7f7f"
colorGray = "#cccccc"
colorGreen = "#11ff00"

all_Drive4 <- read.csv('../../../data/TT1/preprocessed/All/TT1_Drive_4_1s_1s.csv')
all_Drive4$Subject <- as.factor(all_Drive4$Subject)
all_Drive4$logPerspiration <- log(all_Drive4$Perspiration)

persons         = c("01", "02", "03", "04", "05", "06", "07", "09", "12", "13", "15", "16", "17", "18", "22", "24", "29", "30", "31", "32", "41")
starting_points = c( 62 ,  83 ,  70 ,  69 ,  71 ,  61 ,  69 ,  74 ,  65 ,  61 ,  64 ,  63 ,  74 ,  65 ,  79 ,  64 ,  69 ,  69 ,  67 ,  65 ,  62 )
peak_points     = c( 65 ,  86 ,  73 ,  73 ,  73 ,  64 ,  73 ,  79 ,  69 ,  64 ,  68 ,  67 ,  77 ,  68 ,  82 ,  67 ,  72 ,  72 ,  71 ,  68 ,  64 )
ending_points   = c( 67 ,  89 ,  75 ,  76 ,  75 ,  67 ,  75 ,  82 ,  74 ,  66 ,  73 ,  70 ,  80 ,  71 ,  85 ,  70 ,  76 ,  75 ,  75 ,  70 ,  66 )

# Driving time
driving_times = vector(mode="list", length = length(persons))
names(driving_times) <- persons

complete_times = vector(mode="list", length = length(persons))
names(complete_times) <- persons

data_baseline = vector(mode="list", length=length(persons))
pp_baseline = vector(mode="list", length=length(persons))

names(data_baseline) <- persons
names(pp_baseline) <- persons

# Number of peaks

TRACKIN_DURATION = 10
DRIVE_MODE = 4

idx <- 1
plt_AllAcc <- vector(mode="list", length=length(persons)) 
names(plt_AllAcc) <- persons

COLOR_ACC = "#02A3C8"
COLOR_PP = "#F28E8E"
COLOR_BRAKE = "#888888"

y1 <- list(
  tickfont = list(color = COLOR_ACC),
  title="Degree",
  range=c(0, 100)
)
y2 <- list(
  tickfont = list(color = COLOR_PP),
  overlaying = "y",
  side = "right",
  title = "Log Perspiration",
  showgrid = FALSE,
  range=c(min(all_Drive4$ppLogNormalized), max(all_Drive4$ppLogNormalized))
)
  
for (p in persons) {
  pData <- all_Drive4[all_Drive4$Subject==as.integer(p) | all_Drive4$Subject==p,]
  
  # Baseline
  data_baseline[[p]] <- read.csv(str_interp("../../../data/TT1/preprocessed/T0${person}/T0${person}_Drive_1.csv", list(person=p)))
  # Compute the mean
  p_pp_nr <- data_baseline[[p]]$Perspiration
  p_pp_nr <- p_pp_nr[!is.na(p_pp_nr)]
  pp_baseline[[p]] <- log(mean(p_pp_nr))
  
  # Incident
  driving_times[[p]] <- max(pData$Time)
  complete_times[[p]] <- ifelse(starting_points[idx] + TRACKIN_DURATION > driving_times[[p]], driving_times[[p]], starting_points[idx] + TRACKIN_DURATION)
  
  incident_starting_time <- starting_points[idx]
  incident_ending_time <- ending_points[idx]
  complete_time <- complete_times[[p]]
  
  from_time <- ifelse(incident_starting_time - 30 >= 0, incident_starting_time - 30, 0)
  to_time <- complete_time
  
    
  pDataBefore <- pData[pData$Time < incident_starting_time & pData$Time >= from_time,]
  pDataAfter <- pData[pData$Time >= (incident_starting_time + 3) & pData$Time <= to_time,]
  
  ppMeanBefore <- mean(pDataBefore$ppLogNormalized)
  ppMeanAfter <- mean(pDataAfter$ppLogNormalized)
  
  dir.create(file.path('../figures/drive/', paste0('Drive_', DRIVE_MODE)), showWarnings = FALSE)
  fname <- str_interp('../figures/drive/Drive_${drive}/P${person}.svg', list(drive=DRIVE_MODE, person=p)) 
  
  pData <- pData[pData$Time >= from_time & pData$Time <= to_time,]
  plot_Acc <- plot_ly(pData, x = ~Time) %>%
              add_trace(name="Acceleration", y = ~Acceleration, type = 'scatter', mode = 'lines', line=list(width=1.5, color=COLOR_ACC)) %>% 
              add_trace(name="Brake", y = ~Braking, type = 'scatter', mode = 'lines', line=list(width=1.5, color=COLOR_BRAKE)) %>%
              add_trace(name="PP", y = ~ppLogNormalized, type = 'scatter', mode = 'lines', line=list(width=1.5, color=COLOR_PP), yaxis = "y2") %>% 
              add_segments(x = min(pData$Time), xend = max(pData$Time), y = ppMeanBefore, yend = ppMeanBefore, 
                           yaxis = "y2", name="Mean PP (Before Incident)",
                           line=list(color=COLOR_PP, dash = 'dot')) %>%
              add_segments(x = min(pData$Time), xend = max(pData$Time), y = ppMeanAfter, yend = ppMeanAfter, 
                           yaxis = "y2", name="Mean PP (After Incident)",
                           line=list(color="darkred", dash = 'dot')) %>%
              # add_segments(x = min(pData$Time) - 0.1, xend = max(pData$Time), y = pp_baseline[[p]], yend = pp_baseline[[p]], 
              #              yaxis = "y2", name="Baseline PP (from Drive 1)",
              #              line=list(color="blue", dash = 'dot')) %>%
    
              layout(
                title=paste("Subject #", p), 
                xaxis=list(title="Time [s]"), 
                yaxis=y1 , 
                yaxis2=y2, 
                margin = list(l = 50, r = 50, b = 50, t = 50, pad = 4),
                shapes = list(
                  list(type = "rect", fillcolor = "red", 
                       line = list(color = "red"), opacity = 0.3,
                      x0 = starting_points[idx], x1 = ending_points[idx], xref = "x",
                      y0 = 0, y1 = 100, yref = "y"),
                  list(type = "rect", fillcolor = "pink", 
                       line = list(color = "pink"), opacity = 0.3,
                      x0 = starting_points[idx], x1 = complete_times[[p]], xref = "x",
                      y0 = 0, y1 = 100, yref = "y")
                ),
                legend = list(x = min(pData$Time), y = 60)
              )
  
  # orca(plot_PP, fname)
  idx <- idx + 1
  plt_AllAcc[[p]] <- plot_Acc
}

htmltools::tagList(plt_AllAcc)

idx <- 1
behavioralColumns <- c("Subject", 
                       "Brake_u", 
                       "Brake_std", 
                       "PP_before",
                       "PP_u",  
                       "PP_std",
                       "PP_dev")
behavioralMatrix <- matrix(nrow=length(persons), ncol = length(behavioralColumns))

# Careful about Subject 09
# selected_persons <- persons[persons != "09"]

for (p in persons) {
  pData <- all_Drive4[all_Drive4$Subject==as.integer(p) | all_Drive4$Subject==p,]
  
  starting_time <- starting_points[idx]
  ending_time <- ending_points[idx]
  complete_time <- complete_times[[p]] 
  
  from_time <- ifelse(starting_time - 30 >= 0, starting_time - 30, 0)
  to_time <- complete_time
    
  dfBefore <- pData[pData$Time < starting_points[idx] & pData$Time >= from_time,]
  dfAfter <- pData[pData$Time >= starting_points[idx] + 3 & pData$Time <= to_time,]
  
  # diffSpeed <- mean(dfAfter$Speed) - mean(dfBefore$Speed)
  brakeMean <- mean(dfAfter$Braking)
  brakeStd <- mean(dfAfter$Braking)
  
  ppMean <- mean(dfAfter$ppLogNormalized)
  ppBefore <- mean(dfBefore$ppLogNormalized)
  ppStd <- sd(dfAfter$ppLogNormalized)
  
  mid_avg <- (pp_baseline[[p]] + mean(dfBefore$ppLogNormalized)) / 2
    
  diffPP <- mean(dfAfter$ppLogNormalized) - mean(dfBefore$ppLogNormalized)
  
  behavioralMatrix[idx, ] <- c(p, 
                               round(brakeMean, digits=5), 
                               round(brakeStd, digits=5),
                               round(ppBefore, digits=5),
                               round(ppMean, digits = 5),
                               round(ppStd, digits=5),
                               round(diffPP, digits=5))
  idx <- idx + 1
}

# behavioralMatrix

behavioralDf <- as.data.frame(behavioralMatrix)
names(behavioralDf) <- behavioralColumns

behavioralDf

NA

clusteringDf <- behavioralDf
clusteringDf$Subject <- NULL
# clusteringDf$PP_dev_norm <- as.numeric(clusteringDf$PP_dev) / as.numeric(clusteringDf$PP_u)
# clusteringDf$PP_std_norm <- as.numeric(clusteringDf$PP_std) / as.numeric(clusteringDf$PP_u)
clusteringDf$Brake_u <- NULL
clusteringDf$Brake_std <- NULL
clusteringDf$PP_before <- NULL
clusteringDf$PP_u <- NULL
clusteringDf$PP_std <- NULL
# clusteringDf$PP_dev <- NULL

rownames(clusteringDf) <- paste0("#", persons)

for (col in names(clusteringDf)) {
  clusteringDf[,col] <- as.numeric(as.character(clusteringDf[, col]))
  clusteringDf[,col] <- scale(clusteringDf[,col])
}
clusteringDf

library(dendextend)

NUMBER_OF_CLUSTERS = 3

color_darkpink = "#e75480"
CLUSTER_BRANCH_COLORS <- c("red", color_darkpink, "blue")[1:NUMBER_OF_CLUSTERS]
CLUSTER_LABEL_COLORS <- c("red", color_darkpink, "blue")[1:NUMBER_OF_CLUSTERS]

behavioralMatrixClustering <- as.matrix(clusteringDf)
rownames(behavioralMatrixClustering) <- paste0("#", persons)
distMatrix <- dist(behavioralMatrixClustering)
hresults <- distMatrix %>% hclust

hc <- hresults %>% 
      as.dendrogram %>%
      set("nodes_cex", NUMBER_OF_CLUSTERS) %>%
      set("labels_col", value = CLUSTER_LABEL_COLORS, k=NUMBER_OF_CLUSTERS) %>%
      # set("leaves_pch", 19) %>%
      # set("leaves_col", value = c("gray"), k=NUMBER_OF_CLUSTERS) %>%    
      set("branches_k_color", value=CLUSTER_BRANCH_COLORS, k=NUMBER_OF_CLUSTERS)

plot(hc)
legend("topright", 
     legend = c("Exceptional Chicken" , "Normal Chicken" , "No-change/Cool"), 
     col = c("red", "pink" , "blue"), 
     pch = c(20,20,20), bty = "n",  pt.cex = 1.5, cex = 0.8 , 
     text.col = "black", horiz = FALSE, inset = c(0, 0.1))

dfActivity <- all_Drive4[all_Drive4$Time==60,] %>% select(c("Subject", "Activity"))
getActivityName <- function(x) {
  if(x == 1) return("Normal")
  if(x == 2) return("Cognitive")
  if(x == 3) return("Motoric")
}
dfActivity$ActivityName <- sapply(dfActivity$Activity, getActivityName)
dfActivity$Subject <- as.factor(dfActivity$Subject)
rownames(dfActivity) <- NULL
dfActivity %>% select(c("Subject", "ActivityName"))

library(cluster)
fit <- kmeans(clusteringDf, 2)
clusplot(clusteringDf, fit$cluster, color=TRUE, shade=TRUE,
   labels=2, lines=0)

silhouette_score <- function(k){
  km <- kmeans(clusteringDf, centers = k, nstart=25)
  ss <- silhouette(km$cluster, dist(clusteringDf))
  mean(ss[, 3])
}
k <- 2:10
avg_sil <- sapply(k, silhouette_score)
plot(k, type='b', avg_sil, xlab='Number of clusters', ylab='Average Silhouette Scores', frame=FALSE)

# plot(clusteringDf$PP_Dev)
# ggplot(clusteringDf, aes(PP_Mean)) + 
#   geom_density(alpha = 0.3)

# plot(clusteringDf$PP_Dev)
# ggplot(clusteringDf, aes(PP_Dev)) + 
#   geom_density(alpha = 0.3)

ML Model

df_Drive3 <- read.csv('../data/output/Drive_3/corr_Prev_15s_Next_5s.csv')
df_Drive3$Subject <- as.factor(df_Drive3$Subject)

CLUSTER_THRESHOLD = 4
clusters <- cutree(hresults, h=CLUSTER_THRESHOLD)
getClusterName <- function(s) {
  sID <- str_replace(s, "Subject ", "")
  if (str_sub(sID[1], 1,1) != "#") {
    sID <- paste0("#", sID)
  }
  return(ifelse(clusters[sID] == 1, 1, 0))
}

colClass <- sapply(persons, getClusterName)
df_Drive3$clsPP <- as.factor(colClass)

# install.packages("randomForest")
# install.packages("MLmetrics")
# install.packages("caret")
library(randomForest)
library(ROCR)
library(caret)

# Train
trainAndTestModel <- function(pData) {
  pSelected <- pData
  pSelected$Subject <- NULL
  # pSelected$`Steering..std..` <- NULL
  
  # View(pSelected)
  
  nChicken = nrow(pSelected[pSelected$clsPP == 1,])
  nNormal = nrow(pSelected[pSelected$clsPP == 0,])
  
  print(paste("Chicken =", nChicken))
  print(paste("Normal =", nNormal))
  
  # Split dataset
  set.seed(43)
  n_folds = 4
  
  folds <- createFolds(factor(pSelected$clsPP), k = n_folds, list = FALSE)
  # folds <- cut(seq(1,nrow(pSelected)), breaks=n_folds,labels=FALSE)
  pSelected$fold <- folds
  
  kf_acc <- double(n_folds)
  kf_prec <- double(n_folds)
  kf_recall <- double(n_folds)
  kf_f1 <- double(n_folds)
  kf_npv <- double(n_folds)
  
  for (i in 1:n_folds) {
    # train_ind <- sample(seq_len(nrow(pSelected)), size = smp_size)
    # train <- pSelected[train_ind, ]
    # test <- pSelected[-train_ind, ]
    
    # test_ind <- which(folds==i, arr.ind=T)
    train <- pSelected[pSelected$fold != i, ] %>% select(-fold)
    test <- pSelected[pSelected$fold == i, ] %>% select(-fold)
    
    print(paste("Train =", nrow(train), "Pos =", nrow(train[train$clsPP == 1,]), "Neg =", nrow(train[train$clsPP == 0,])))
    print(paste("Test =", nrow(test), "Pos =", nrow(test[test$clsPP == 1,]), "Neg =", nrow(test[test$clsPP == 0,])))
    
    # Model Train
    model <- randomForest(clsPP~., data = train, importance = TRUE, ntree=10)
    print(importance(model))
    # print(model)
    
    # Test
    testX <- select(test, -clsPP)
    predY <- predict(model, testX)
    testY <- test$clsPP
    
    #print(levels(predY))
    #print(levels(testY))
    
    # Evaluate
    # print(table(predY, testY))
    
    kf_acc[i] <- mean(predY==testY)
    kf_recall[i] <- sensitivity(predY, testY)
    kf_prec[i] <- posPredValue(predY, testY, positive = 1)
    kf_f1[i] <- (2 * kf_recall[i] * kf_prec[i]) / (kf_recall[i] + kf_prec[i])
    kf_npv[i] <- negPredValue(predY, testY, positive = 1)
    
    # print(paste("Perf:", kf_acc[i], kf_recall[i], kf_prec[i], kf_f1[i], kf_npv[i]))
  }
  
  # XGB
  param <- list(objective       = "binary:logistic", 
               booster          = "gbtree",
               eval_metric      = "auc",
               eta              = 0.1,
               max_depth        = 5,
               gamma            =0.8,
               min_child_weight = 3,
               subsample        = 1,
               colsample_bytree = 0.5,
               stratified       = F
  )
  pSelected <- pSelected %>% mutate(clsPP=ifelse(clsPP==1, 1, 0))
  
  # AUC
  # aucs = c()
  # xgb_m = xgb.cv(   params               = param,
  #                   data = as.matrix(pSelected %>% select(-clsPP)) ,
  #                   label =  pSelected$clsPP,
  #                   nrounds             = 500,
  #                   verbose             = F,
  #                   prediction          = T,
  #                   maximize            = T,
  #                   nfold = n_folds,
  #                   metrics  = "auc",
  #                   early_stopping_rounds = 100,
  #                   scale_pos_weight = 1)
  # aucs =  c(aucs,as.numeric(xgb_m$evaluation_log[xgb_m$best_iteration,"test_auc_mean"]))
  
  # Get average performance
  acc <- mean(kf_acc)
  prec <- mean(kf_prec)
  recall <- mean(kf_recall)
  f1 <- mean(kf_f1)
  npv <- mean(kf_npv)
  # auc <- mean(aucs)
 
  # Return 
  rtn <- list(
    accuracy=acc,
    recall=recall,
    precision=prec,
    f1=f1,
    npv=npv
    # auc=auc
  )
  
  print(rtn)
}

trainAndTestModel(df_Drive3)

library(caret)
library(VGAM)

NUM_FEATURES = 8

pSelected <- df_Drive3
pSelected$Subject <- NULL

fit <- vglm(clsPP~., family=multinomial, data=pSelected)
# summarize the fit
summary(fit)

probabilities <- predict(fit, pSelected[,1:NUM_FEATURES], type="response")
predictions <- apply(probabilities, 1, which.max)
predictions[which(predictions=="High")] <- levels(pSelected$clsPP)[1]
predictions[which(predictions=="Low")] <- levels(pSelected$clsPP)[2]
# summarize accuracy
table(predictions, pSelected$clsPP)

# x <- pSelected[,1:NUM_FEATURES]
# y <- pSelected[,NUM_FEATURES + 1]
# # fit model
# fit <- plsda(x, y, probMethod="Bayes")
# # summarize the fit
# summary(fit)
# # make predictions
# predictions <- predict(fit, pSelected[,1:NUM_FEATURES])
# # summarize accuracy
# table(predictions, pSelected$clsPP)

library(caret)

pSelected <- df_Drive3
pSelected$Subject <- NULL
pSelected$clsPP <- as.numeric(colClass)

# define training control
train_control <- trainControl(method = "cv", number = 4)

# train the model on training set
model <- train(clsPP ~ .,
               data = pSelected,
               trControl = train_control,
               method = "glm",
               family=binomial())

# print cv scores
summary(model)

predTrain = predict(model, newdata=pSelected, type="raw")
table(pSelected$clsPP, predTrain > 0.5)

formula = clsPP ~ .
modellm <- glm(formula, data=pSelected)
plot(modellm)
summary(model)